Fractal Threshold Behavior in Vacuum Gravitational Collapse 
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We present the numerical evidence for fractal threshold behavior in the five dimensional vacuum 
Einstein equations satisfying the cohomogeneity-two triaxial Bianchi type-IX ansatz. In other words, 
we show that a flip of the wings of a butterfly may influence the process of the black hole formation. 
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Introduction. — Critical phenomena in gravitational 
collapse are an interdisciplinary field that was initiated 
by the remarkable work of Choptuik [1]. It concerns 
the study of the basin boundary between two generic 
states — space-times with and without black holes. This 
bistable behavior shares many properties with the phase 
transition in statistical mechanics. It is a well-known fact 
that, in general, basin boundaries can be either smooth 
or fractal, and there is no reason to believe that gravity 
is special in this context. However, in more than one 
hundred papers that were devoted to critical phenomena 
in gravitational collapse so far [2, 3] , the basin boundary 
between dispersion and collapse to a black hole is always 
smooth, and there is no indication of chaos [3| . The single 
counterexample [6] is restricted to the solution with many 
unstable modes and as such it does not give a chance to 
observe fractal threshold behavior directly in the initial 
value problem. 

The theory of chaos in infinite dimensional systems is 
not fully developed yet, but this field is of great interest 
and growing rapidly. In this Letter we provide an exam- 
ple of chaos in reduced Einstein equations that constitute 
a system of partial differential equations. We clarify 
the hypothesis of Bizoii et al. ^ and present fractal 
threshold behavior that can be interpreted as a chaotic 
scattering in a critical surface — the surface in the phase 
space of initial data that separate collapse to a black hole 
from dispersion. This surface is smooth, but it contains 
three copies of the critical solution, and the basin sets 
of these copies have fractal boundaries. In this sense, 
the chaotic behavior presented here is double critical and 
can be directly seen in the dynamical evolution. In our 
Letter we estimate the fractal dimension of the basin 
boundaries. The fractal dimension is a diffeomorphism 
invariant indicator of chaos in general relativity [8| . 

Setting. — We consider vacuum gravitational collapse 
in a very simple setting. Namely, in order to evade 
Birkhoff's theorem and save radial symmetry, we take 
five dimensional vacuum Einstein equations and reduce 
the number of degrees of freedom by the BCS ansatz 
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where A, 5, B, and C are functions of time t and radius 
r. One-forms (jfc are standard left-invariant one-forms on 
SU{2), 
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where < 6* < tt, < (/) < 27r, and < i/" < 47r are Euler 
angles. The angular part of the metric ([T]) corresponds 
to the spatial part of the metric in Bianchi IX cosmology. 
SU{2) is dipheomorphic to S^; therefore B and C func- 
tions can be interpreted as squashing parameters of S^. 
B ^ C = correspond to the round and B = C ^ 
to the deformed with one squashing parameter (the 
so-called biaxial case). In this Letter we consider the 
triaxial case where B, C are independent functions — 
two dynamical degrees of freedom. The vacuum Einstein 
equations derived for the metric H]) (see f^) possess a 
discrete symmetry. It corresponds to the freedom 
of permutations of coefficients of one-forms ak in the 
angular part of the metric ([T]). These permutations are 
generated by the following transpositions: 
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where the transposition swaps the coefficients of erf 
and crj in ([1]). Biaxial configurations correspond to the 
fixed points of these transpositions: {B,B), {B,~B/2), 
and {B,—2B). All biaxial solutions exist in three geo- 
metrically equivalent copies. _ 

Critical phenomena. — The study of the critical 
phenomena in the model restricted to biaxial symmetry 
revealed type-II gravitational collapse with the critical 
solution being discretely self-similar (DSS). It was also 
shown Q that in triaxial symmetry the phenomenology 
of critical behavior remains the same. The Z3 symmetry 
implies that in the triaxial case the critical solution 
(hereafter denoted as DSSi — because it is discretely self- 
similar and has one unstable mode) takes one of three 
equivalent forms. Let us denote them as x[^\ Xj^'' 
and ^3^'. Each of these copies acts as a codimension- 
one attractor in the phase space of solutions. Hence, 
the codimension-one critical surface Merit contains three 
basin sets Mi {i — 1,2,3). It was shown by Bizon 
et al. [7] that the basin boundaries of these sets are 
given by the stable manifolds of the discretely self-similar 
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codimension-two attractor DSS2- Moreover, the authors 
suggested that the basin boundaries could be fractal. 

In this Letter we continue the work that was initi- 
ated in Q and investigate the structure of the basin 
boundaries of codimension-one attractors X^'^'K In other 
words, we study numerically the Cauchy problem for 
two-parameter families of initial data. We fine-tune one 
of the parameters to the critical surface A4crit and the 
remaining parameter to the basin boundary (this means 
fine-tuning to the codimension-two attractor DSS2)- 

One of the problems in carrying on such studies 
numerically comes from the fact that in order to cancel 
two growing modes of DSS2 in general two-parameter 
initial data, one is forced to use very high numerical 
precision. It leads, on the software level, to the necessity 
of multiprecision packages. These packages slow down 
the numerical evolution, and the calculations become 
computationally infcasible. The possible solution is to 
find special two-parameter initial data that are conve- 
nient in the study of the basin boundaries in Aicrit- We 
propose the following procedure for finding such initial 
data. 

Let us consider general initial data that are fine-tuned 
to DSS2- The solution generated by such initial data 
approach DSS2 and, since it is always fine-tuned with 
a finite precision, repels from DSS2- However, one 
can freeze the evolution at the moment the solution is 
"nearly" DSS2- Using this snapshot one can generate 
another two-parameter family of initial data. 

In particular, we consider time-symmetric initial data 
that were studied in 

BiO,r)=pfir), C(0,r)-ai?(0,r), (4) 

where /(r) = 100 cxp[— 20(r — 0.1)^] is the generalized 
Gaussian (for details, see Q)- In order to fine-tune to 
DSS2 we set parameters 

a ^ 0.1411036683285, , , 

V = 0.09524484187150217214296769, ^ ' 

and evolve the solution up to t = 1.038. Next, we 
introduce two new parameters k, i and rescale C{t,r) = 
KC{t, r), B{t, r) — LB{t, r). In this way we create a new, 
two-parameter (k, l) family of initial data. 

Hereafter we restrict our analysis only to this partic- 
ular family. The critical surface Adcrit corresponds to 
a curve l*{k) in the parameter space. We present this 
curve in Fig. [1] 

In order to characterize evolution of initial data given 
by l*{k) we define below a map /i : k e M ^ a G 
{1,— 2, 0}. Generic solutions starting from l*{k) 

approach one of the codimension-one attractors x[^\ 
X2^\ X^^^"^ and recover biaxial symmetry C = aB in one 
of three equivalent forms: a = 1, a = — i, a — —2, 
respectively. Moreover, there are points that separate 
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FIG. 1: The critical surface Merit presented in the phase 
space of parameters k and l. All initial conditions defined by 
(.*(«) lead to the solutions that approach the critical solution 
DSSi or attractors of higher codimension {DSS2 or possibly 
others), i- ~ 1-* leads to supercritical solutions (with a black 
hole) and t ~ t* gives subcritical solutions (dispersion). The 
error bars are lower than 10~^. 

basins of attraction A4i. For a practical convenience 
we denote them as a ^ 0, but we stress that this is 
not related to the scaling C — aB. The set of all such 
points forms the basin boundary Q. It was shown in Q 
that solutions defined by Q approach the DSS2 solution 
or possibly higher codimension attractors. The map 
h : K ^ a describes asymptotic behavior of solutions, 
but in finite precision numerical calculations it describes 
intermediate dynamics. Such calculations show that if 
KA < Kb and h{KA) 7^ h{KB), and h{KA)h{KB) ^ 0, then 
there exists kq such that ka < kq < kb and h{KQ) = 0. 
Again, in practical calculations we mark by /i(k) — 
all points with an undetermined basin set. For these 
points the double precision of the parameter l* is not 
sufficient to cancel by a bisection the A^crit-transversal 
growing mode of DSS2- This mode dominates over the 
A^crit-tangential growing mode of DSS2 and the DSSi 
behavior is not observed. In contrast to this, for /i(k) ^ 
the yVfcrit-tangential growing mode of DSS2 takes over 
and "pushes" solutions along Merit- Such solutions 
approach one of the copies of DSSi. Finally, since the 
bisection has a finite precision, they are also repelled from 
Aicrit along the unstable mode of DSSi. This behavior 
was described in detail in [3]. 

To study the geometry of Merit, especially the ge- 
ometry of basin sets A4i, we determined the map h 
in the numerical experiment. We have collected two 
samples corresponding to two scales in the basin bound- 
ary. Sample 1 contains 15001 points and corresponds to 
K G S*! = [—0.5,1]. Sample 2 contains 3201 points and 
corresponds to k ^ S2 — [0.07,0.078]. All test points 
were uniformly distributed in respective intervals. 

It is clearly seen in Fig. [5] that the basin boundaries 



FIG. 2: Sample 1. The type of the DSSi solution as a function 
of K (plotted as a histogram). h{n) = 1, — —1/2, h{K) — 



—2 correspond to X[^\ Xj^', X!^''' , respectively, and h{K) — 
corresponds to points with an undetermined basin. 



-(1) 



(1) 



have a complex structure Magnifying a part of the 
boundary reveals a new structure presented in Fig. [31 
Therefore, one should expect that an intersection of the 
basin boundary Q with tested sets Si , S2 has nontrivial 
fractal dimension < dim{Sir)Q) < 1, where i = 1,2. In 



order to estimate the fractal dimension we follow [llLll2| 
and calculate the uncertainty dimension [3] (from now 
on dim denotes the uncertainty dimension). 

Let 5' be a one-dimensional set in one-dimensional 
parameter phase space (we have one free parameter k). 
The probability that any two random points ka, kb 
separated by a distance e belong to different basins 
h{KA) ^ Hkb) scales as P(e) - gi-rf™(snQ)^ Testing 
many pairs of points for several values of e one can fit 
P{e) to the power law and estimate dim{S n Q). The 
wider the range of e we consider (assuming e is small 
enough), the higher precision we obtain. 

The fit to the power law for Sample 1 presented in Fig. 
Hreveals that dim{Si HQ) = 0.722 ± 0.006, therefore the 
basin boundary is fractal. However, the quantity esti- 
mated above is an "averaged" fractal dimension. In order 
to see this, let us consider a set S{k) = [k — 0.2, k + 0.2] 
and dim(S{K) n Q) for k such that S{k) C Si. It follows 
from Fig.[5]that the uncertainty dimension dim[S{K)C]Q) 
varies with k, so the investigated structure seems to be a 
multifractal rather than a monofractal. Such hypothesis 
is supported by the fractal dimension of Sample 2, namely 
dim{S2r\Q) = 0.680±0.006 (see the linear fit in Fig.©. It 
follows from the definition of the uncertainty dimension 




FIG. 4: Sample 1. P is a probability that any two random 
points KA, i^B separated by a distance e belong to different 
basins Ii^ka) 7^ /i(«:_g). The uncertainty dimension was 
estimated to be dim{Si n Q) = 1 - w = 0.722 ± 0.006, where 
w is & slope of the linear fit above. 



that 5*2 C 5i ^ dim{Si n Q) > dim{S2 n Q). The 
estimated values satisfy this inequality as expected. The 
smaller e we take, the more precise bisection we need 
and more quickly the numerical double precision of t* 
is exhausted. Therefore, Sample 2 contains ca. 18% of 
points with an undetermined basin in contrast to only 
ca. 7% in case of Sample 1. Clearly, the finite precision 
of the bisection makes the number of points with an 
undetermined basin scale dependent and may have a 
greater effect on the fractal dimension of Sample 2. The 
error bars shown in Figs. |4] and [6] indicate statistical 
errors 
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The multifractal objects can be more completely char- 




FIG. 3: Sample 2. Part of Fig. [2] magnified. The complex 
structure is seen at a different scale. 



FIG. 5: The uncertainty dimension dim{S{K)r\Q) varies with 
K. The set S{it) — [n — 0.2, k + 0.2] contains 4001 points for 
each value of k. The horizontal dotted line indicates averaged 
uncertainty dimension dim{Si HQ). The gray area represents 
errors. 
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FIG. 6: Sample 2. The uncertainty dimension was estimated 
to be dim{S2 n Q) = 1 — w — 0.680 ± 0.006, where lo is a slope 
of the hnear fit above. 

acterized by the singularity spectrum, but such analysis 
of our system would involve more precise data and an 
enormous computational power. Let us mention that 
the computer resources needed to collect Sample 1 and 
Sample 2 involved over six dozens of processors with the 
integration time measured in weeks. 

The direct analogy to finite dimensional dynamical 
systems suggests that the DSS2 solution plays a role 
of a strange repeller but we do not have sufficient 
dynamical evidence to support such claim. To the best 
of our knowledge, the nature of solutions that drive the 
dynamics in a fractal basin boundary remains unknown 
in the context of partial differential equations A 
construction of the DSS2 solution in ^ is a first example 
that may give some insight into this problem. 

Conclusions. — In summary, we have presented the nu- 
merical results that indicate the presence of fractal basin 
boundaries in the critical surface in five dimensional vac- 
uum gravitational collapse. Moreover, the data support 
the hypothesis that basin boundaries are multifractals. 
We note that all three copies of the critical solution {x[^\ 
X2^^ and ^3^'') are geometrically equivalent. Therefore, 
one should not expect to see chaos in the black hole mass 
scaling of supercritical solutions. The solutions behave 
chaotically in the nearly double-critical regime, and then 
settle down to the Minkowski space-time or the black 
hole. This final state is not sensitive to initial conditions, 
but the way in which it is approached is, indeed, chaotic. 
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